clear
clc
close all
A=-4;
utv=@(t) exp(-10.*t.^2);
g1uv=@(u) u.^2;
g2uv=@(u) A.*u;
  
syms t u;
yu=utv(t);
dyu=diff(yu,t);
gv=dyu-g1uv(utv(t))-g2uv(utv(t));
f1=g1uv+gv;
f2=g2uv/u;
f=f1+f2*u;
f1tu=matlabFunction(f1);
% f2t=matlabFunction(f2);
f2t=A;
ftu=matlabFunction(f);
% gt=matlabFunction(gv);
ut=matlabFunction(yu);


t1=0;
t2=1;
u0=ut(t1);

%  N=[100 1000 10000];
N=[64 128 256 512 1024];
ujn=[];

for i=1:length(N)
    [ujn,err,xx,uex]=odesemi_R_k3th(t1,t2,N(i),u0,ut,f1tu,f2t,ftu);
    %error:L2error,L∞error
    h1=(t2-t1)/N(i);
    
    L2err(i)=norm(err,2).*h1.^0.5;
    Lwuerr(i)=norm(err,Inf);
    figure(1); semilogy(xx,err); hold on;
    figure(2);plot(xx,err);
    set(gca,'FontSize',10)
    xlabel('x','FontSize',12)
    ylabel('|u(x_i,1)-u_i^1|','FontSize',12)
end
%ORDER
L2order=log(L2err(1:end-1)./L2err(2:end))./log(N(1:end-1)./N(2:end));
Lwuorder=log(Lwuerr(1:end-1)./Lwuerr(2:end))./log(N(1:end-1)./N(2:end));
fprintf('explicit 4_Order  Runge-Kutta,\n')
fprintf('ode\n')
fprintf('du/dt=f(t,u)+g(t) t∈(t_1,t_2)\n')
L_G_output(N,L2err,Lwuerr,L2order,Lwuorder);



function [ujn,err,xt,uex]=odesemi_R_k3th(t1,t2,T,u0,ut,f1tu,f2t,ftu)

tau=(t2-t1)/T;
xt=linspace(t1,t2,T+1);
ujn(1)=u0;

for iw=1:T
    ti=xt(iw);
    yi=ujn(iw);
    k0=ftu(ti,yi);
    k1=(f1tu(ti+2/3*tau,yi+2/3*tau*k0)+f2t*yi)...
        /(1-2/3*tau*f2t);
    k2=(f1tu(ti+2/3*tau,yi+2/3*tau*k1)+f2t*(yi+1/6*tau*k1))...
        /(1-1/2*tau*f2t);
    k3=(f1tu(ti,yi-1/2*tau*k1+1/2*tau*k2)+f2t*(yi-3/2*tau*k1+1/2*tau*k2))...
        /(1-tau*f2t);
    ujn(iw+1)=ujn(iw)+tau/8*(3*k1+3*k2+2*k3);
    

end

uex=ut(xt);
err=abs(ujn-uex);


end

